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ABSTRACT 


The newly launched X-ray satellite, eRosita, has recently revealed two gigantic bubbles extending to ~ 80° above and below the 
Galactic center. The morphology of these “eRosita bubbles” bears a remarkable resemblance to the Fermi bubbles previously dis- 
covered by the Fermi Gamma-ray Space Telescope and its counterpart, the microwave haze. The physical origin of these striking 
structures has been intensely debated; however, because of their symmetry about the Galactic center, they likely originate from 
some energetic outbursts from the Galactic center in the past. Here we propose a theoretical model in which the eRosita bubbles, 
Fermi bubbles, and the microwave haze could be simultaneously explained by a single event of jet activity from the central super- 
massive black hole a few million years ago. Using numerical simulations, we show that this model could successfully reproduce 
the morphology and multi-wavelength spectra of the observed bubbles and haze, which allows us to derive critical constraints on 
the energetics and timescales of the outburst. This study serves as an important step forward in our understanding of the past 
Galactic center activity of our Milky Way Galaxy, and may bring valuable insights into the broader picture of supermassive black 


hole-galaxy co-evolution in the context of galaxy formation. 


Introduction 


Recent data taken by the eRosita satellite has revealed striking images of two giant X-ray bubbles extending ~ 80° (which 
corresponds to ~ 15 kpc assuming a distance to the Galactic Center (GC) of 8.5kpc) above and below the GC!. Despite the 
larger extent, the morphology of the “eRosita bubbles” is remarkably similar to the Fermi bubbles, the two gamma-ray bubbles 
detected by the Fermi Gamma-ray Space Telescope in 20107. It was also shown previously that the gamma-ray bubbles have 


counterparts in the microwave band, known as the microwave/WMAP/Planck haze***, as well as polarized lobes in radio”. The 


edges of the Fermi bubbles at low latitudes also align with earlier detections of X-ray shells in the 1.5 keV band by ROSAT®’. 
Because of their morphological similarities, enormous physical sizes, and symmetry about the GC, these fascinating structures 
most likely originate from the same event of powerful energy outburst from the GC sometime in the past. Understanding their 


physical origin would therefore provide valuable information about the history of our Milky Way Galaxy. 


Before the detection of the eRosita bubbles, the physical origin of the Fermi bubbles and microwave haze had been hotly 
debated®. Key questions to address include whether the gamma-ray emission is hadronic or leptonic, and whether the triggering 
mechanism is associated with a nuclear starburst, or with activity of the central supermassive black hole (SMBH). Because of 
the proximity to the GC, the ample spatially-resolved, multi-wavelength observational data provide very stringent constraints on 
the proposed theoretical models. One strong constraint comes from the hard spectrum of the gamma-ray bubbles with spectral 
index of —2.1°, which means that the gamma-ray emission has to be generated by cosmic-ray (CR) protons transported by 
starburst or active galactic nucleus (AGN) winds that do not easily cool (i.e., the “hadronic wind models")!!! by CR electrons 


quickly transported by AGN jets before they cool (i.e., the “leptonic jet models")!7-! 


, or by CRs which are accelerated in situ 
by shocks or turbulence near the shock fronts (i.e., the “in situ acceleration models")!®!8_ Previous studies have found that 
purely hadronic models are disfavored because the predicted microwave emission from secondary particles produced via pion 
decay is not enough to account for the observed haze emission’. In addition, the upper limits in TeV gamma-rays obtained by 
HAWC have ruled out hadronic models with a power law extending to PeV energies!’. On the other hand, both the leptonic 
jet models and in situ acceleration models remain promising mechanisms to explain the gamma-ray bubbles and microwave 
haze!” 17-18. Tn this work, we show that the new eRosita data provides crucial information that allows us to put additional 


constraints on these two scenarios, and that the combination of the gamma-ray, X-ray, and microwave images and spectra 


strongly suggest that past jet activity of the GC black hole is the likely culprit. 


Simulations of Past Jet Activity of Sgr A* 


We performed three-dimensional hydrodynamic simulations of energy release at the GC in the form of bipolar jets perpendicular 
to the Galactic plane. The simulations self-consistently model the evolution of cosmic rays (CRs) injected with the jets at the 
GC, including their dynamical interactions with the thermal gas within the Galactic halo, and energy losses of CR electrons 
due to synchrotron radiation and inverse-Compton (IC) scattering as they travel within the Galactic magnetic and radiation 
fields (see Methods for the simulation details). Modeling of the CRs together with the thermal gas allows us to compute the 
thermal radiation from the gas and the non-thermal emission generated by the CRs self-consistently. For the simulation results 
presented in this work, the jet activity of the central SMBH is assumed to start 2.6 Myr ago and last for 0.1 Myr (see Methods 
for discussion on the chosen parameters). Because of the large pressure contrast with respect to the ambient gas, the jet material 
expanded into a pair of “cocoons” or “bubbles" above and below the GC, similar to radio bubbles observed in galaxy clusters. 
At the present day, the cocoons have grown and reached a height of ~ 7.5 kpc from the Galactic plane (Figure 1). The CR 


electrons within the cocoons that were transported from the GC interact with the interstellar radiation field (SRF) and shine in 
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the gamma-ray band as the observed Fermi bubbles extending to Galactic latitudes of |b| ~ 50° — 55° (Figure 2). The same 
energy injection from the black hole and subsequent cocoon expansion pushed the gas within the Galactic halo away from 
the GC with supersonic speeds, forming an outward propagating shock. At the shock front, the compression of gas caused an 
increase in the local gas density, producing enhanced thermal Bremsstrahlung emission in the X-ray band manifested as the 


eRosita bubbles. 


As can be seen from Figure 2, the morphology and sizes of both the eRosita bubbles and Fermi bubbles are successfully 
reproduced by the leptonic jet model. Importantly, unique features including their brightness distributions, sharp edges, and 
smooth surfaces match observational constraints well (see Figure 3 and Methods for detailed discussion). For the Fermi bubbles, 
reproducing all the above characteristics simultaneously has proved to be a challenging task for many theoretical models®, 
because they require very specific conditions to be met during the bubble formation. For instance, in the leptonic scenario, the 
intensity of the ISRF decays with increasing distance from the GC. Consequently, the CR energy density needs to be somewhat 
enhanced toward higher Galactic latitudes. This enhancement could be caused by initial adiabatic compression when the AGN 
jets were active!®!5, The sharp edges of the gamma-ray bubbles require suppression of CR diffusion across the bubble surface, 
which could be a result of anisotropic CR diffusion along magnetic field lines that are wrapped around the bubble surface due 
to the magnetic draping effect!*. Furthermore, the smooth surface of the Fermi bubbles indicate that there are no vigorous 
hydrodynamic instabilities on the scale of the bubble sizes. This could be explained by insufficient time for the instabilities to 


grow!?, or by suppression of instabilities owing to viscosity”. 


The X-ray emission predicted by the leptonic jet model shows very good agreement with the observed eRosita bubbles 
as well, not only in terms of the extension of the X-ray bubbles, but also in terms of the X-ray surface brightness variations. 
Figure 4 shows the comparisons between the simulated and observed X-ray surface brightness profiles at three horizontal cuts, 
|b| = 40°, 50°, 60°, as well as the latitude-averaged profiles. Overall, the predicted amplitudes of the brightness variations and 
the locations of the X-ray bubble edges are largely consistent with the observational data. In our simulation, the forward shock 
compresses the gas into a thick shell with width of ~ 2.5 kpc (Figure 1), consistent with estimates from simple geometrical 
models!. Due to the projection of the compressed gas shell enclosing the gamma-ray bubbles with low gas densities, the 
modeled X-ray profiles are limb-brightened at the Galactic longitude of |/| ~ 50°, similar to the data. At lower latitudes 
(|b| = 40° and 50°), there is some emission at |/| ~ 15° — 20° due to gas that is ejected within the AGN jets being compressed 
near the contact discontinuity. This feature may be related to the two peaks at |/| ~ 20° for |b| = 40° seen in the observed 
profile for the north X-ray bubble. Despite the large scatter seen in the data, the simulated latitude-averaged profile of the X-ray 
emission as a function of Galactic longitude approximately falls between the observed values corresponding to the emission 
from north and south bubbles. It is difficult to establish an exact match to the data though, because the observed X-ray sky in 
this region has many complex features. In particular, the northeast part of the eRosita bubble is coincident with a prominent 
structure called North Polar Spur (NPS). The NPS was discovered several decades ago, but its origin remains elusive and is a 


subject of ongoing debate. Due to the spatial correlation with the Loop I feature in the radio band, some proposed that the NPS 
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is associated with a local superbubble in the Solar neighborhood*!;*”; other evidence supports scenarios related to past GC 
activity”? ?5. The discovery of the south eRosita bubble has further strengthened the GC hypothesis. Regardless of its origin, 
the NPS generates a substantial degree of asymmetry in the X-ray sky, complicating the interpretations of the true emission 
from a GC event. While these considerations prohibit a firm conclusion, our simulated X-ray profiles suggest that an energy 
outburst from the GC could contribute to the NPS emission, at least partially (i.e., the peak of emission at / ~ —20°). This is in 
line with the conclusion drawn from a recent study of the radio/optical polarization data near the NPS*°, which suggests that 


the NPS could be a superposition of both local and GC structures (see Methods for more detailed discussion). 


Additional constraints on the formation of these gigantic bubbles come from the broad-band spectra in the same region in 


k>-4, GeV gamma rays by Fermi”, to TeV gamma rays by HAWC!”. We compiled 


the sky from microwaves by WMAP and Planc 
the available observational constraints and our simulated spectra in Figure 5. In agreement with previous studies’, we find 
that the leptonic model can simultaneously reproduce the GeV gamma-ray spectrum of the Fermi bubbles and the microwave 
spectrum of the WMAP/Planck haze. In the model, the primary CR electrons injected by the AGN jets are quickly transported 
to large distances before they cool due to synchrotron and IC losses. Therefore, the CR electrons could maintain their hard 
spectrum (spectral index of -2.1) during propagation, which is necessary for producing the hard gamma-ray and microwave 
spectra seen in the data. In addition, there appears to be a high-energy cutoff in the Fermi data at ~ 110 GeV, and so far no 
detection in the TeV range has been made. This high-energy cutoff is also a generic feature predicted by the leptonic jet model. 
In the early phases of the evolution, the CR electrons suffer great synchrotron and IC losses due to strong magnetic field and 
radiation field near the GC, generating an exponential cutoff at high energies in their spectrum. After the jets leave the GC, the 
dynamical timescale of the jets becomes shorter than the synchrotron and IC cooling times, and therefore the CR spectral shape 


is essentially frozen during the later expansion!>. Overall, the broad-band spectrum from the bubbles can be well explained by 


a single population of CR electrons that were transported from the GC by jets from the central SMBH. 


The good agreement between the simulation predictions and the images and spectra of the Fermi bubbles, eRosita bubbles, 
and microwave haze, suggests that past jet activity from the GC black hole could be the common origin of these fascinating 
structures in our Milky Way Galaxy. The jet parameters adopted in the simulations are summarized in Table 1. Although these 
parameters may not be unique (see Methods for discussion), they could still provide us with valuable insights into the formation 
processes. We find that the new X-ray data from eRosita allows better constraints on the duration of the jets, as the duration 
controls the separation between the forward shock (edges of the X-ray bubbles) and the contact discontinuity (edges of the 
gamma-ray bubbles) in the model. Our simulations suggest that the forward shock at the present day has reached a vertical 
height of ~ 11 kpc away from the Galactic plane and a width of ~ 14 kpc, i.e., an oblate ellipsoid rather than a sphere (Figure 
1). Because the sizes of the eRosita bubbles are comparable to the distance between the Sun and the GC (~ 8.5 kpc), the top 
and bottom of the eRosita bubbles in fact correspond to the near side of the shock front which is only ~ 6 kpc away from the 
Sun roughly in the vertical direction, rather than the top of the shock front near the rotational axis of the Galaxy. In the shock 


downstream, the electron number density is enhanced to ne ~ 1073 cm~? and the gas is heated to T ~ 10° K. At the present day, 
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the forward shock is moving at a speed of ~ 2000 km s~! (Mach number M ~ 10) in the vertical direction and ~ 1300 km s7! 
(M ~ 6) in the lateral direction at a height of 5 kpc away from the Galactic plane. Inside the contact discontinuity, the electron 
number density is ne ~ 1075 — 1074 cm™? and gas temperature is T ~ 10° — 10’ K (though the temperature distribution within 
the contact discontinuity is uncertain due to parameter degeneracies; see Methods). According to our model, the central SMBH 
was active ~ 2.6 Myr ago, injecting a pair of bipolar jets in mostly kinetic forms for a duration of ~ 0.1 Myr. After taking into 
account uncertainties in the initial conditions (Methods), the Sgr A* was estimated to be accreting at ~ 0.1 — 1 the Eddington 


rate during the active phase, corresponding to a consumption of ~ 10° — 104 solar masses within ~ 0.1 Myr. 


Our model of the Fermi/eRosita bubbles makes a specific prediction that could be tested with future observations. The region 


3 and 


between the contact discontinuity and the forward shock is compressed to electron number densities of a few ~ 107°cm7 
shock heated to a few 10’ K, which are the conditions comparable to those typically observed in the intracluster medium of cool 
core galaxy clusters. Unlike in most of the Milky Way halo, this region could produce iron Kæ emission from hydrogen-like 
(6.7 keV) or helium-like (6.9 keV) iron ions that have peak collisional ionization equilibrium fractions at a few 10’ K and 
108 K, respectively”’. We note that whether such a line is detectable depends on the fraction of the fast hot gas. Since the 
temperature of the gas is uncertain in our model, that fraction may be too small for existing instruments to detect. However, it 
may be detectable with Athena as its effective area is almost two orders of magnitude larger than Chandra’s at 6.9 keV7*. Given 
the expected post-shock velocity of ~ 2000 km s~! (see lower right panel in Figure 1), the total Doppler line broadening is 
expected to be about 45 eV at 6.7 keV line energy, but the exact value will depend on the details of the velocity projection. A 
shift of this magnitude should be detectable with the Athena X-ray Integral Field Unit that has a planned spectral resolution of 
2.5 eV up to line energies of 7 keV and more than sufficient spatial resolution to probe this region’. Thus, a coherent outflow 


originating from the region between the shock and contact discontinuity could be a testable model prediction, especially if the 


data quality is sufficient to constrain not only the line width but also its skewness or profile. 


Discussion 


As mentioned in the Introduction, both the leptonic jet model and the model of in situ acceleration by shocks or turbulence 
are promising mechanisms to explain the gamma-ray emission of the Fermi bubbles, but new insights could be obtained by 
adding the constraints from the eRosita data. One primary difference between the leptonic jet model and the shock acceleration 
model is the location of the forward shock driven by the outburst. In the leptonic jet model, the surface of the eRosita and 
Fermi bubbles corresponds to the forward shock and contact discontinuity, respectively. On the other hand, if the Fermi bubbles 
are produced by shock accelerated CRs as proposed in some of the in situ acceleration models*”, the more extended eRosita 
bubbles would be left unexplained. In this regard, the leptonic jet model may be favored as both the eRosita bubbles and 
the Fermi bubbles could be simultaneously accounted for by a single event. In situ acceleration models in which the Fermi 
bubbles are produced by CRs accelerated by turbulence generated within the bubbles!* remain a possibility. In this scenario, 


the forward shock driven by the outflow could account for the eRosita bubbles, and the turbulence in the shock downstream 
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could potentially stochastically accelerate CRs that are responsible for the gamma-ray bubbles. However, note that the surface 
of the Fermi bubbles is very smooth, suggesting that hydrodynamic instabilities at the bubble surface are suppressed, and hence 
the amount of turbulence generated may be limited. Also, in this scenario the turbulence is expected to be volume-filling in the 
shock downstream and of increasing strength toward the shock front!*, and therefore the sharp edges of the Fermi bubbles 
cannot be naturally explained. More detailed simulations are needed to test whether the stochastic acceleration models could 
satisfy all the observational constraints. 

Our simulation suggests that a single past jet activity ~ 2.6 Myr ago from the GC black hole could be the common origin 
for the eRosita bubbles, the Fermi bubbles, and the microwave haze. The required Eddington ratio is estimated to be ~ 0.1 — 1 
during the brief active accretion phase of ~ 0.1 Myr. Though such high activity of Sgr A* may appear somewhat surprising 
given its extreme quiescence at the present day, multiple lines of evidence have pointed to much elevated past activity of Sgr 
A* over the past ~ 10 Myr®!. In particular, recent studies have found enhanced ionization levels in the Magellanic Stream, 
which could be explained if Sgr A* went through a phase of Seyfert-like flaring activity a few Myr ago’. The energetics 
and timescales obtained by our simulations are in good agreement with such a Seyfert-flare scenario. The epoch of outburst 
predicted by our model is also in line with several independent constraints from kinematic studies of the halo gas using X-ray 


and UV absorption lines*?-*4 


, which estimated that the outburst occurred ~ 2 — 8 Myr ago (though there remain uncertainties 
in the interpretations; see Methods for detailed discussion). 

The rich multi-wavelength observational data as well as detailed theoretical modeling have provided valuable information 
about the past GC activity. For example, our model suggests that the radiation and magnetic fields are likely suppressed in 
the inner kpc region near the GC at the time when the jets were launched (see Methods). Such suppression could perhaps 
be related to an earlier injection that produced the GC chimney?’ and bipolar radio bubbles*° recently observed close to the 
Galactic plane. Alternatively, the suppression could be associated with the formation of the nuclear star cluster ~ 6 Myr ago>’, 
where stellar winds or supernova explosions of the massive young stars could have helped to evacuate the gas and create a 
GC environment less hostile to the cooling CR electrons within the jets. Such a scenario is also broadly consistent with the 
emerging picture of SMBH-galaxy co-evolution in that AGN activity often accompanies star formation activity in the galaxy*®. 
Some open questions remain to be addressed, e.g., whether CRs accelerated at the forward shock could generate observable 
non-thermal emission in the leptonic jet scenario, whether the thermal structure of the Galactic halo probed by X-ray absorption 
line studies?’ could be explained, how the jet-induced outflows entrain cold gas and how it may be related to the high velocity 
clouds observed in the Galactic halo*’. Future investigations will further reveal the impact of this energetic feedback on the 
evolution history of our Milky Way Galaxy, and how this event fits in the broader picture of SMBH-galaxy co-evolution in the 


universe. 


6/23 


Methods 


Simulation setup and parameters 

We carried out three-dimensional hydrodynamic simulations of bipolar jets emanating from the GC and perpendicular to the 
Galactic plane including relevant CR physics using the FLASH code!*:*°-4!, We utilized the CRSPEC code to self-consistently 
model the evolution of CR spectrum due to synchrotron and IC cooling, while accounting for the dynamical interaction between 
the CRs and the thermal gas!>. The modeling of the CR spectral evolution is crucial, as it allows us to simultaneously compute 
the non-thermal radiation from the CRs and the Bremsstrahlung emission from the thermal gas. The injected CR electrons are 
assumed to follow an initial power-law spectrum ranging from 10 GeV to 10 TeV, with a spectral index of -2.1. In our previous 
work!*, we demonstrated that the magnetic field within the Fermi bubbles needs to be amplified to values comparable to the 
ambient field at the present day in order to reproduce the microwave haze emission. Therefore, we do not simulate the magnetic 
field directly, but adopt the default magnetic field distribution in GALPROP” for the computation of the microwave haze, 


i.e., 


B| = Boexp(—z/zo)exp(—R/Ro), where R is the projected radius to the Galaxy’s rotational axis, Bo is the average field 
strength at the GC, and zg and Ro are the characteristic scales in the vertical and radial directions, respectively. We adopt zo = 2 
kpc, Ro = 10 kpc, and Bo = 50uG as motivated by observational constraints**. As will be discussed later, the high-energy 
cutoff in the observed gamma-ray spectrum requires suppression of radiation and/or magnetic field strengths near the GC at the 
time of jet injection. Therefore, for computing the synchrotron losses in the simulations, we have used a factor of three lower 
normalization of magnetic field strength compared to the value adopted in our earlier work!>. For IC scattering, we utilize the 
ISRF model in GALPROP and compute the CR electron energy losses and gamma-ray emissivity including the Klein-Nishina 
effect. To generate the X-ray image, the X-ray emissivity in the energy range 0.6-1.0 keV is calculated using the MEKAL 
model implemented in XSPEC. The initial condition for the halo gas is the same as in previous leptonic jet models!% 13, which 
assumes an isothermal halo with T = 2 x 10° K in hydrostatic equilibrium within a fixed Galactic potential. The normalization 
of the gas density profile is chosen to match the observed profile derived from X-ray absorption line studies*>. Other simulation 


13-15 


setups are essentially identical to our previous works , though the new eRosita data motivated some parameter adjustments. 


We summarize the parameters in Table 1 and discuss key differences as follows. 


In previous studies of the leptonic jet model!” !° 


, it was found that there are parameter degeneracies when the morphology 
of the gamma-ray bubbles was the primary constraint. The new eRosita data provides another key constraint on the separation 
between the forward shock and the contact discontinuity in the model. Compared to the parameters adopted in our previous 
simulations!>, we find that the duration of the jets has to be reduced from tiet = 0.3 Myr to tje = 0.1 Myr. Changing the jet 
duration alone would produce Fermi bubbles that are too oblate. In order to maintain the right axial ratio of the Fermi bubbles, 
we then need to decrease the total (thermal plus CR) injected energy density to one third of the original value, i.e., from 
ej + ejer = 4.09 x 107° to 1.36 x 107° erg cm~>. By changing these two parameters, the morphology of both the X-ray and 


gamma-ray bubbles can be reproduced. Note, however, that the overall dynamics and bubble shapes are controlled by the total 


injected energy density (thus the total pressure within the bubbles), and hence the individual values of ej and ejcr are degenerate 
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in the current model. Due to this parameter degeneracy, the predicted gas temperature within the Fermi bubbles should be 
considered uncertain as it depends on the amount of injected thermal energy assumed in the model. Additional observational 


33,34, 46,47 


constraints such as X-ray, UV, and radio emission/absorption lines may be used to break the parameter degeneracy. 


We will investigate this issue in future studies. 


The simulations are analyzed using the software y°', and the simulated gas and CR distributions at t = 2.6 Myr are shown 
in Figure 1. At this time, the forward shock driven by the AGN outburst has propagated to a vertical height of ~ 11 kpc anda 
radius of ~ 7 kpc from the rotational axis of the Galaxy. The shock front is currently expanding with a speed of ~ 1000 — 2000 
km s~!, with nearly vertical flows near the axis and progressively larger opening angles away from the axis. In the shock 
downstream, the gas is compressed to ne ~ 107°? cm~? and heated to T ~ 108 K, producing the X-ray bubbles as observed. 
Within the contact discontinuity is the underdense cavity filled with thermal gas and CRs, though their relative contributions are 
uncertain due to the parameter degeneracy discussed above. The CR distribution is edge-brightened, which is necessary for 
reproducing the nearly flat intensity of the Fermi bubbles!>. Similar to our previous findings, we find that only a small fraction 
(fe ~ 107°) of the simulated CRs needs to be CR electrons in order to match the observed gamma-ray and microwave spectra. 
Note that the value of fẹ is larger than that found in our previous studies because of the reduced amount of CR energy injected 


in the current simulations. 


For completeness, we show the comparisons between the predicted and observed profiles for the microwave haze and 
gamma-ray bubbles in Figure 2. As can be seen, the leptonic jet model not only can simultaneously reproduce the spectra of the 
gamma-ray bubbles and the microwave haze (Figure 3), but also the key features in their spatial distributions. The predicted 
synchrotron emission follows a centrally peaked profile, primarily because of the exponential decay of the Galactic magnetic 
field from the GC. As for the gamma-ray profiles, the nearly flat intensity profiles and sharp edges of the observed bubbles are 
reproduced, which requires very specific three-dimensional distributions of CRs as well as suppression of CR diffusion across 
the bubble surface (see discussion in the main text). Note, though, that due to large levels of noise and asymmetries in the data, 
it is impossible to provide a perfect fit to the data. Therefore, we do not aim to complicate the models by introducing additional 
free parameters to improve the fits as it could result in overinterpretation of the data. However, should future observations 
be characterized by lower noise, such improvements could be warranted. For instance, the mismatch at / ~ 20° between the 
predicted and observed gamma-ray profiles may be related to the bending of the observed Fermi bubbles toward the west and 
that could be accounted for in the model by invoking a slight jet tilt as considered in our earlier work!. As a side note, we point 
out that in our previous analysis based on magnetohydrodynamic simulations'*, the CRs needed to be replenished by a second 
jet injection in order to match the microwave profile. That was because the initial magnetic field strength in the simulation was 
set to its present-day value, and hence the strong magnetic pressure near the GC pushed the gas and CRs away and caused a 
depression in the microwave profile at r < 15°. However, our current analysis, which is based on hydrodynamic simulations 
that do not model the effects of a dynamically important initial magnetic field, shows that a single jet could also yield consistent 


48 


results with the observed haze profile. Since the second jet scenariot was no longer supported by the updated Fermi data’, this 
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would imply that a single jet plus an initially suppressed magnetic field is the more likely scenario. As we will show later, the 


constraints from the high-energy cutoff in the gamma-ray spectrum also point to the same conclusion. 


Comparisons with other observational constraints 

Our model predicts the existence of a thick shell of ultrahot gas (T ~ 108 K) moving with high velocities (~ 1000 — 2000 km 
s7!) within the Galactic halo. Is this consistent with existing X-ray observations of the Galactic halo and kinematics of the halo 
gas as traced by X-ray/UV absorption lines from background quasars? While the outflow velocities and timescales of outburst 


32-34 other studies have inferred more mild outflow 


predicted in our model are consistent with some of the observational studies 
velocities of ~ 200 — 300 km s7! 17-49-50, The apparent discrepancies could be explained by taking into account a number of 
factors: (1) the predicted hot gas in the range of T > 108 K would not contribute substantial emission in the X-ray band probed 
by current data. In addition, for strong shocks, it is likely that the ions are preferentially heated to higher temperatures than 
the radiating electrons. Therefore, outflow velocities inferred by measuring the temperature contrast (where the temperature 
is obtained by fitting the X-ray spectrum emitted by the electrons) between the shocked gas and the ambient medium may 
underestimate the true velocities*’; (2) the UV absorption lines probe the kinematics of the cooler T ~ 104 — 10° K gas, and 
thus using them to infer the velocity of the hot gas would depend on how exactly the cooler gas is formed/entrained in the 
hot flow, which remains a major unresolved question itself. Furthermore, in an entrainment scenario, the velocities of the 
cold clouds are typically smaller than those of the hot gas due to imperfect momentum transfer. The inferred velocities of 
~ 1000 — 1300 km s~! from UV absorption lines?™’4 would thus imply even higher velocities for the hot flow; and (3) there 
remain large uncertainties in the interpretations of the observational data, e.g., assumptions about the geometry of the outflows 
and injection patterns (continuous vs. instantaneous injections), the asymmetric emission measure of the Galactic halo*!, and 
confusion due to foreground/background projections (e.g., the NPS). Given all the above considerations, we leave the detailed 


comparisons with these observational data to dedicated future work. 


Initial GC environment and particle acceleration 

As discussed in our previous study!>, the latitude-independent high-energy cutoff in the observed gamma-ray spectrum at 
~ 110 GeV implies a very uniform distribution for the maximum CR energy, Emax ~ 300 GeV. The maximum energy of CR 
electrons today is set by the fast synchrotron and IC cooling near the GC at early times, followed by adiabatic cooling due to 
the bubble expansion. These considerations allow us to infer conditions within ~ kpc from the GC in the early phase of the 
evolution. Because of the parameter modifications mentioned above, the expansion is somewhat slower, and the age of the 
bubbles becomes somewhat longer than our previous estimates, changing from 1.2 Myr to 2.6 Myr. This leads to longer times 
for the CR electrons to cool adiabatically, changing Emax by a factor of 10 from the simulation time t ~ 0.4 Myr to t = 2.6 Myr. 
This implies that the maximum CR energy is ~ 3 TeV at the time when the jets leave the inner kpc region. According to the 
constraints we obtained in the previous study! (see their Figure 6), it implies that either the initial jet velocity needs to be larger, 


or the initial radiation plus magnetic field energy density near the GC needs to be suppressed when the jets were first launched. 
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Assuming the initial jet velocity is unchanged (otherwise it would modify other properties such as bubble morphology), the 
radiation plus magnetic field energy densities inferred would need to be suppressed to Utot = UB + Urag ~ 1.9 x 1071? erg cm~?. 
Assuming the energy density of the ISRF at the GC is urag ~ 1.54 x 107!? erg cm~3, the same as the GALPROP model at the 
present day after taking into account the Klein-Nishina effect, then the magnetic field strength at the GC at the time of injection 


would need to be ~ 3 uG. Equivalently, if one would like to estimate characteristic CR cooling times, one could define an 


effective magnetic field strength such that utot = |Befr|*/ (87). Our constraints would imply |Befe| ~ 74G. 


One could ask whether the inferred GC conditions are consistent with CR electrons being accelerated up to ~ 3 TeV 
by checking two criteria. First, assuming the CRs within the jets are accelerated by diffuse shock acceleration (DSA), the 
maximum CR electron energy can be estimated by balancing the DSA acceleration rate and the IC plus synchrotron cooling 


rate>2, 


B -1/2 
Ee max ~ 6 X 10*E!/2 ($4) TeV, 


where € ~ 0.1 is the DSA acceleration efficiency. One can see that, given the suppressed radiation and magnetic field strengths, 
the acceleration of CR electrons is not strongly limited by the synchrotron and IC cooling. Of course, the magnetic field 
strength within the jets is uncertain and could be higher, especially close to the black hole*>. But the above estimate suggests 
that the condition of Ee max > 3 TeV could be satisfied for a wide range of magnetic field strength. The second criterion is the 
constraint on acceleration timescales. The maximum rate of particle acceleration by a strong parallel shock can be written in 


the following form**>°, 


dE 2 B _ 
7 7 15x10 7 (5) Vi.GeVs 1, 


where Z is the charge, and V,, is the shock velocity in cm s~!. For Vey = 0.025c, B ~ 3uG, and Z = 1, dE /dt ~ 8.44 x 1077 
GeV s~!. One could then estimate the required acceleration timescale face = Emax /(dE /dt) to accelerate CRs to Emax = 3 
TeV, which is much less than the dynamical time of the jets of fayn ~ (1 kpc) /(0.025c) ~ 0.13 Myr. In other words, there is 
sufficient time for CR electrons to be accelerated to ~ 3 TeV. Therefore, overall we find that DSA is a plausible mechanism to 


accelerate CR electrons to the energy required in the model. 


Relativistic magnetic reconnection (RMR) is also a possible primary particle acceleration mechanism. The properties of 
RMR depend strongly on the magnetization parameter o = B? / (4rpc?) (here, p is the rest mass particle energy density), with 
the relativistic regime defined by o > 1. It is found that as o increases, the energetic particle spectral index p (assuming the 
particle energy spectrum is a power law distribution, E~?) decreases, and can drop below p = 2. This raises the interesting 
possibility that the observed E~? spectrum is actually an aged spectrum, which would allow the jet to be older and/or the 
particle transport slower than what is inferred from DSA. However, the hardening of the spectrum is accompanied by a decrease 


in the maximum particle energy achievable (as it must be): Ymax ~ [(o6 + 1)(2 — p)/(p—1)]!/2-?). Thus, to achieve 3 TeV 
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(Ymax ~ 6 x 10°), either o must be very large or p must be very close to 2. For example, if p = 1.5, o ~ 2.4 x 103, while if 
p = 1.9, o could be as small as 42. But with p so close to 2, the constraints on aging are essentially the same as for DSA. Thus, 


we regard the constraints on jet age and particle transport as robust. 


Uncertainties and limitations of the current model 
The power of the two jets obtained in our model is 6.32 x 10+ erg s~!, which is close to the Eddington luminosity for Sgr A*. 
Note, however, that the estimated jet power in our model is directly proportional to the unknown initial central gas density 
assumed in the simulations!*. Although we have used the observational constraints from X-ray absorption line studies* to 
inform the gas profile of the Galactic halo, the true value near the GC at the time of injection remains highly uncertain. As 
discussed in the main text, if the gas near the GC were evacuated due to a prior energy injection or the formation of the nuclear 
star cluster, it is conceivable that the central density is lower than that assumed in our model. Therefore, the energetics of the 
SMBH accretion event obtained in our model shall be considered as an overestimate, and the Eddington ratio is likely to be in 
the range of ~ 0.1 — 1 given the above considerations. 

Our simulation setup of the Galactic halo is relatively simple and is strictly symmetric about the GC. In reality, the Galactic 
halo is much more complex. The most prominent structure is the NPS, which is very bright in the north-eastern X-ray sky. 
As mentioned in the main text, it remains uncertain whether the NPS is associated with a local superbubble or a GC event. 


21,22,26 whereas 


Generally speaking, studies based stellar polarization and extinction tend to support the local bubble scenario 
analyses on the basis of X-ray data tend to suggest a GC origin?*. In order to reconcile the apparent discrepancies, it has 
been proposed that the X-ray emission of the NPS may be a superposition of both the local and GC structures”°. Our simulation 
also supports this conclusion. It has also been reported that the emission measure of the Galactic halo is asymmetric about the 
Galactic plane>!. If the outflow from the GC is expanding into this asymmetric Galactic medium, the brighter emission of the 
NPS could also be accounted for”. Another limitation of our simulation is that we have only modeled the hot component of the 
halo gas and neglected the colder, multi-phase Galactic disk. Therefore, our simulation results cannot be directly compared 
to observed structures close to the Galactic plane, including the X-ray shells observed by ROSAT®, and the more recently 
discovered HI outflows*®, X-ray chimney” and bipolar radio bubbles*©. Since all these structures may provide important 
clues to the formation of the gamma-ray/X-ray bubbles, we will extend our model to include the Galactic disk component and 
investigate the detailed jet-disk interactions in future work. 

In order to inflate the nearly symmetric eRosita and Fermi bubbles, our model requires bipolar jets in the direction 
perpendicular to the Galactic plane. Generally speaking, the directions of SMBH jets are determined by the orientation of 
accretion disks and/or the black hole spins on much smaller scales and do not need to align with the rotational axis of the host 
galaxies. Future studies would be required to see whether this condition could be relaxed by considering tilted jets interacting 
with dense, multiphase interstellar medium within the Galactic disk>®. Alternatively, such alignment may indeed be expected in 
low-mass, disky galaxies like the Milky Way, as shown by recent simulations>’. In this case, the rotational axis of the black 


hole accretion disk may align with that of the host galaxy as it is fed by gas with high angular momentum, and the black hole 
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spin could efficiently align with the accretion disk due to the Bardeen-Petterson effect because of the relatively small mass of 
the black hole compared to its accretion disk. 

Finally, we note that the computational costs of our three-dimensional, CR hydrodynamic simulations forbid us to do a full 
parameter scan of all the possible combinations of the six jet parameters. Therefore, the jet parameters we found in this study 
may not be unique, even though the parameter space should be very limited given the stringent constraints from all the available 
observational data!*. Nevertheless, our current simulation serves as a proof of concept that the gigantic bubbles within our 


Galaxy could plausibly originate from past jet activity of the GC black hole. 


Data Availability 
Simulation data that were used to generate the emissivity profiles are available in the Supplementary Information. Source data 


associated with other main figures of the manuscript are available from the corresponding author upon reasonable request. 


Code Availability 
The simulations were performed using the code FLASH, publicly available at https://flash.rochester.edu/site/flashcode/, with 
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Figure 1. Simulated gas and CR properties. Slices of the simulated (a) gas density, (b) temperature, (c) CR energy density, 
and (d) absolute magnitudes of the velocity field at t = 2.6 Myr. Streamlines show the directions of the outflows driven by the 
jet injection from the GC. 
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Figure 2. Mock all-sky maps centered at the GC. The morphology of the simulated gamma-ray (dark purple; 1-2 GeV) and 
X-ray (blue-yellow; 0.6-1 keV) bubbles resembles the Fermi bubbles and the eRosita bubbles, respectively. Our simulations 
suggest that these fascinating structures in the Milky Way Galaxy can be generated by a powerful outburst from the GC 
supermassive black hole about 2.6 million years ago. 
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Figure 3. Gamma-ray and microwave profiles. (a) Intensity profiles of the microwave haze at 23 GHz as a function of 
radius from the GC. The blue curve shows the predicted synchrotron emission from CR electrons (CRe) from the leptonic jet 
model, and blue shaded region shows the observed limits of the WMAP haze’. (b) Gamma-ray intensity profiles (N is the 
number of gamma-ray photons per unit area per solid angle per unit time) as a function of Galactic longitudes averaged over the 
latitude region of 40° < |b| < 50°. Curves show the predicted IC emission, and shaded regions show the observed limits’. 
Colors correspond to different energy bins (red: E = 1 — 3 GeV; green: E = 3 — 10 GeV; blue: E = 10 — 500 GeV). 


16/23 


|b| = 60° |b] =50° |b] = 40° Averaged 


2x10! 
a a — Model b c d 
2 — Observed north bubble 
ks] — Observed south bubble 
7 
un 
2 10! 4 4 
S 
n 
n 
ov 
i= 
& 
5 6x10° 
È 
a 
w 
© o 
£ 4x10 
5 
a 

3x 10° r r r r r r r r r r r r r r r r r r 7 z 

-100 -50 to) 50 100 -100 -50 to) 50 100 -100 -50 0 50 100 -100 -50 to) 50 100 
Galactic longitude (deg) Galactic longitude (deg) Galactic longitude (deg) Galactic longitude (deg) 


Figure 4. X-ray surface brightness profiles. Profiles of X-ray surface brightness (in units of counts s~! deg~*) as a function 
of Galactic longitudes in the 0.6-1.0-keV range comparing the model predictions (blue) and the observed profile of the eRosita 
bubbles (north and south bubbles shown in orange and red colors, respectively). Panels from left to right correspond to 
horizontal cuts at Galactic latitudes of |b| = 60°,50°, 40°, and the latitude-averaged profiles. The overall amplitudes and 
locations of the X-ray shells are well reproduced by the model, though interpretation is complicated by the asymmetry in the 
data caused by the North Polar Spur observed in the northern X-ray sky. 
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Figure 5. Broad-band spectrum for the Fermi bubbles. Multi-wavelength spectrum for the region containing the Fermi 
bubbles ranging from microwave emission observed by WMAP and Planck, GeV gamma rays by Fermi, and TeV gamma-ray 
upper limits by HAWC. For direct comparisons, we computed the simulated spectra within the same extracted regions as the 
observed spectra, namely, |/| < 10°, |b| = 20° — 40° for the gamma-ray bubbles, and |/| < 25°, |b| = 10° — 35° for the 
microwave haze. The black and red curves show the predicted emission of the leptonic jet model, in which the primary CR 
electrons can simultaneously produce the observed gamma-ray emission due to IC scattering and the microwave haze emission 
by synchrotron radiation, while satisfying the HAWC constraint. 
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Table 1. Input and Derived Jet Parameters 


Parameter Description Value Unit 

n Density contrast 0.05 - 

Ne Energy density contrast 0.81 - 

Cicer CR energy density 6.82 x 1071! erg cm~? 
Vjet Jet speed 0.025 c 

Riet Radius of cross-section 0.5 kpc 

fjet Duration of injection 0.1 Myr 

Mej Electron number density 0.1 cm”? 

Pj Thermal gas densiy 1.95 x 107% gem? 
ej Thermal energy density 1.29x 107? erg cm~? 
Pre Kinetic power 3.08 x 10% erg s7! 
Ph Thermal power 7.21 x 10% erg s7! 
Poy CR power 3.82 x 10%! erg s i 
Pg Magnetic power 3.31 x 107! erg s7! 
Pet Total power 3.16 x 10% erg s7! 
Ejer! Total injected energy 1.00 x 1057 erg 


1 The total injected energy by both bipolar jets is 2Ejet- 
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